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Large deviations of Lyapunov exponents 

Abstract. Generic dynamical systems have 'typical' Lyapunov exponents, measuring the sensitivity to 
small perturbations of almost all trajectories. A generic system has also trajectories with exceptional values 
of the exponents, corresponding to urmsually stable or chaotic situations. From a more mathematical point 
of view, large deviations of Lyapunov exponents characterize phase-space topological structures such as 
separatriccs, homoclinic trajectories and degenerate tori. Numerically sampling such large deviations using 
the Lyapunov Weighted Dynamics allows one to pinpoint, for example, stable configurations in celestial 
mechanics or collections of interacting chaotic 'breathers' in nonlinear media. Furthermore, we show that 
this numerical method allows one to compute the topological pressure of extended dynamical systems, a 
central quantity in the Thermodynamic of Trajectories of Ruelle. 



PACS numbers: 05.45.-a, 05.10.Gg, 05.40.-a, 05.45.Jn 



1. Introduction 



In many physical systems, the typical situations are not the most interesting, and one is led to look for rare 
trajectories. For example, of all possible initial conditions of eight planets, only a small number lead to a 
stability comparable to the one of the solar system [1, 2, 3]. Exceptional configurations such as resonances and 
separatrices play an important role in this kind of system [1]. Similarly, a study of the transport properties 
of almost-integrable systems requires the knowledge of chaotic layers that are extremely thin, because they 
are the structures responsible for the global diffusion mechanism [4, 5, 6]. Apart from being rare, often these 
relevant trajectories turn out to be unstable, a further source of difficulty. For example, unstable soliton and 
breather modes may be important factors in the transport of energy of Bose-Einstein condensates, and even 
macromolecules [7, 8] . Similarly, the phenomenon of intcrmittency, which has a profound effect in turbulent 
systems, is believed to be generated by localized spatial structures that appear and disappear in the flow [9] . 

Each physical situation will have one or more quantities whose large deviations are relevant. For example, 
for a rare 'rogue' wave, a quantity that is certainly interesting is its energy or its height, for a traffic problem 
the actual flow, and for a planetary orbit its eccentricity or inclination angle. Large deviations of Lyapunov 
exponents are a particular and important case, as has been recognized for years [10, 11, 12], which still 
attracts lots of attention from the dynamical system community nowadays [13, 14, 15, 16]. 

Large Deviations 

For a dynamical system defined by trajectories of n-dimensional variables x(t), consider an observable 
A(x). We wish to study the trajectories such that, for example, dt A{x.) = tAq. In particular, their total 
probability: 

P{A,,T) = (s(^j^ dtA{^)-TA,)) (1) 

A more convenient object to study is the Laplace transform of P{Ao,t) 

Z(a,r) = (e"^-^") = j dAoc'^^^" P{Ao,t) (2) 

Note that a is to Ag what the inverse temperature is to the energy (density) in statistical mechanics: Z{a, r) 
is nothing but a partition function in the space of trajectories and the large deviation formalism thus allows 
one to extend the static formalism of statistical physics to dynamical observables [10, 18]. It is thus natural 
to define 

fi{a,T) = - logZ(a,r) (3) 

T 

which plays the role of a free energy and 



(A p^-^-^oX 

A„(a,r)=;.'(a,r) = VbxT^ (4) 



which plays the role of an order parameter, distinguishing between various classes of trajectories with 
markedly different values of Ao- 

Lyapunov exponents and their large deviations 

Consider two nearby points Xi{t — 0) and X2(i = 0) and their subsequent trajectories Xi(i) and X2(t) 
in a dynamical system. The distance Ui{t) = xi(t) — X2(t) at long times will typically grow as 

ln|xi(t)-X2(i)| -iAi (5) 

The quantity Ai = Ai is a measure of the sensitivity to the initial conditions. Similarly, we may consider 
three non-colinear nearby points, and the two differences Mi(t) = xi(i)— X2(i) andtt2(t) = X2(i)— X3(t). The 
lengths of ui{t) and U2{t) will grow as above, but we may ask about the area |ui A U2I of the parallelogram 
defined by them. This defines a new exponent: 

Injui AU2I ~tA2 =t(Ai + A2) (6) 
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In general, for generic points and p vectors {ui, Up}, the p-volume grows as: 

In |ui A U2 A ... A Upl = ^ In I detju^.u,,}] - tAp = t(Ai + A2 + ... + Xp) (7) 

The Aa are the (finite-time) Lyapunov exponents whose large deviations will be the focus of this article. 
This means that we shall consider intervals of time r, and the role of Aq above will be played by, for example, 
Aq = Xi. We can then classify trajectories according to the value of their first Lyapunov exponent A° and 
construct the probability distribution: 

P{X,)^{6{X,~X1)) (8) 

and the corresponding moment- and cumulant-generating functions 

Z(a,r) = (e™^i) /i(a, r) = log Z(a, r) (9) 

and their generalizations to higher Lyapunov exponents. 

At this stage one may ask what is special about large deviations of Lyapunov exponents, that this 
should merit a separate article. There are a number of reasons why we think this is so. To begin with, it 
is clear that Lyapunov exponents are a measure of chaoticity, and rare trajectories having unusually low or 
unusually high values will represent havens of stability (important for, e.g., planetary systems), or chaotic 
bursts. We shall see several examples of this below. 

The sum of the positive Lyapunov exponent is also related (via Pcsin's theorem [17, 18]) to the 
Kolmogorov-Sinai entropy: a space-time quantity 'counting' the number of typical trajectories, a useful 
measure of the geometry of the system. 

Another important measure is the dimension of attractors. These are related to the cumulative Lyapunov 
exponents A^. Consider the value of p such that Ap > and Ap+i < 0. The dimension of the attractor has 
to be intermediate between p and 1, so that a volume element of the attractor adverted by the dynamics 
neither overfiows nor contracts on it: the precise value is encoded in the Kaplan- Yorke formula [17, 18]. 

An even more specific property of Lyapunov large deviations is that it allows us to pinpoint important 
topological features in the dynamics. For instance, the family 

Zp(a = l,r)-(e-^'')-(e-^^:-i^") (10) 

is related to normally hyperbolic invariant manyfold (NHIM) with p unstable directions [19]. These are 
manifolds which are unvariant under the dynamics and whose normal directions have the structures of saddles, 
with exactly p unstable directions. To see this, consider an overdamped Langevin dynamics with very small 
noise. Computing Zi amounts to selecting trajectories that emerge from a saddle point with exactly one 
unstable direction (see figure 1), the changes in distribution due to the acceleration or deceleration along 
the trajectory is exactly compensated by the factor Ai. All degrees of freedom but one lie at the bottom of 
energy wells and Zi has stabilized the unstable manifold emanating from a NHIM with one unstable direction. 
Similarly, computing Z2, where we bias with Ai + A2, exactly compensates the area variations incurred when 
leaving from a saddle with two unstable directions, all other degrees of freedom again lying at the bottom 
of energy wells. In general, computing Zp, the bias with Ap, stabilizes the p-dimensional unstable manifold 
emerging from critical points with p unstable directions. What we have described is the basis of the relation 
between supersymmetry and Morse theory ([20, 21]), expressed in the context of stochastic equations. For 
dynamics more general than purely dissipative, the structures stabilized by Zp are not restricted to unstable 
manifold emerging from critical points, as we shall see below for Hamiltonian dynamics. 

Another area where fluctuations of Lyapunov exponents are expected to play a crucial role is that of 
intermittent dynamics. In such cases, the coexistence in phase space of trajectories with different chaoticities 
can be seen as the signature of a first order phase transition. Measuring the free energy /i(a,<) and order 
parameter A(a) in intermittent system is thus an important challenge. To understand why intermittency can 
be related to first order dynamical transitions, consider the following example. Suppose there is a metastable 
structure, with a lifetime T. Within this structure, the largest Lyapunov exponent A takes an average value 
Xms- On the other hand, the average of A over the typical trajectories (the stable attractor) is, say. As. 

X To maintain the introduction as free of mathematical technicalities as possible, we reserve more precise definitions of (finite- 
time) Lyapunov exponents for section 7 
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Figure 1. Lyapunov exponents and topological structures. The volume contractions and expansions 
induced by evolution are exactly the instantaneous growth of the sums of the p largest Lyapunov exponents. 
Top: The weight exp(Ait) stabilizes the unstable manifold emerging from critical points with one unstable 
direction. Bottom: The weight exp[(Ai + X2)t] stabilizes the unstable manifold emerging from critical 
points with two unstable directions. 

We wish to calculate the large deviation function (e"*'^), and to estimate which trajectories contribute (see 
[22]). If a{Xms — As) > 0, the trajectories that belong to the metastable structure will be favored by the bias 
a in the measure. The quantities aX^s and aXs are the biases per unit time attributed to trajectories in 
metastable and stable basins, they have dimensions time~^. This bias is compensated by the time of escape 
from the metastable state T, so that the condition for the trajectories in the metastable state to dominate 
the large deviation computation is: 

ra(A™s-As)>l (11) 

If {Xms — Xs) is extensive, and the metastable state is rather stable (T is large), the transition takes place 
close to a = [2.3]. It will be sharp, and first order. 

Outline 

In what follows we shall present several applications of the Lyapunov Weighted Dynamics (LWD), a 
cloning algorithm that allows to simulate a population of trajectories (or clones) with a biased measure 
P(A, T)e"'*''^ [24]. We will briefly review the numerical method in section 2 while the technical details, the 
underlying formalism and its extension to the case of several Lyapunov exponents are presented last, in 
section 7, to make the article more readable. 

Then, we will show how the LWD can localize tiny regular island in almost completely chaotic systems 
(section 3) and detect the disappearance of the last regular islands steming from the Lagrange points in 
the restricted gravitational three-body problem (section 4). We will then show in section 5 that regular 
Hamiltonian systems develop positive Lyapunov exponents when perturbed by a weak, additive noise and 
that trajectories with atypical Lyapunov exponents are localized on interesting objects such as separatrix. 
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degenerate tori or unstable manifolds. We will then turn to systems with large number of degrees of 
freedom, showing that trajectories with atypically large kth Lyapunov exponent in the Fermi-Pasta-Ulam 
chain typically involve k chaotic breathers (see section 6.1). Last, we will show in section 6.2 that the 
LWD allows one to compute the dynamical free energy fJ.{a) and order parameter A(a) in one-dimensional 
coupled-map lattices. 



2. Lyapunov Weigthed Dynamics - presentation of the algorithm 

For sake of clarity, let us consider a given one-dimensional stochastic dynamical system 

X = f[x{t)] (12) 

where f{x) is a complicated function that contains some noise terms. Let us call e a parameter controlling 
the noise intensity so that e = is the deterministic limit of the system. To compute Lyapunov exponents, 
one also need to consider the linearized dynamics 

u^f'[x{t)]u (13) 

From the solution of (13), one can then compute the largest Lyapunov exponent 

The LWD is a population Monte Carlo algorithm whose purpose is to sample trajectories according to 
a modified measure so that a trajectory x{t) with a finite-time Lyapunov exponent Xi{t) is sampled with a 
probability 

^a[a;(i)] = ^T^F[x(t)]e"*^i« with Z{a,t) = [ V[x{t)]P[x{t)]e'''^^^'^ (15) 
Z(a,t) J 

where P[x{t)] is the probability to observe the trajectory x{t) with the unbiased dynamics (12). The 
algorithm was introduced in [24] and briefly reviewed in [25] (see also [26] where it was used to localize 
atypical trajectories in convex billiards). It is presented in full details in section 7.3 after a quite extensive 
discussion of the underlying formalism at the beginning of section 7. Nevertheless, before we present several 
applications of the algorithm in the sections 3 to 6, we briefly present the method below and try to give 
hints about why it is efficient and where are its pitfalls. 

The principle of the algorithm is as follows. One simulates Nc copies {xi{t),Ui{t)) of the system, also 
called "clones", with the unbiased dynamics (12) and (13). Then, at every time step, one computes for each 
clone i 

\ui{t - dt)\ 

The clone i is then replaced, on average, by Silt)" copies. If the dynamical system were deterministic (e = 0), 
all the Si(i)" copies would then follow the same trajectory and a given clone would yield, after a time t, a 
total number of copies given by: 

t / dt t / dt , , , s 1 1 / ^ 1 

nil nil Min-mr |U,(0)|" ^ ^ 

If one were to pick a clone at random in the population at time t, the probability to choose a trajectory x{t) 
would thus be proportional to e°'^'^'^P[x{t)], where P[a;(t)] is the probability that the dynamics (12) produces 
the trajectory x{t). This thus results in the biased measure (15). Of course, simulating many times the same 
trajectory is not enhancing the quality of the sampling and will not reveal any rare trajectories that brute- 
force sampling would miss. This is where the stochasticity of f[x{t)] comes into play: when e ^ 0, all copies 
of a given clone made at a time t' will follow different trajectories for t > t'. Since the cloning is larger for 
trajectories that carry a large weight, their vicinity is well sampled by the algorithm. Conversely, trajectories 
with small weights rapidly die out and no computing power is spent on their simulation. In practice, there 
are some differences — mostly due to numerical efficiency — between the LWD and the algorithm presented 
above (for instance, the number of clones Nc is kept constant in LWD), which are detailed in section 7.3. 
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In various parts of this paper, noise is not part of the original problem, and it is added only to let the 
different clones diffuse independently, to enhance the sampling of the algorithm. In those cases, the level 
of noise e has to be taken as small as possible and one should ensure that the e — > limit is not singular. 
Finally, e can be taken to zero at the end of the run to check that the atypical trajectories that have been 
discovered are indeed solution of the equations of motion of the original system. 

When the system conserves energy, it may be desired to search for trajectories within a specified energy 
shell - and the same can be said of other constants of motion, if present. It is then useful to consider 
energy-conserving noise. This may be implemented in a number of ways, such as rotating the velocity vector 
randomly, or, more generally proposing a small random change of the coordinates and projecting them 
back onto the energy shell (see Appendix A for a momentum- and energy-conserving noise in Hamiltonian 
systems) . 

The values of a and e play roles analogous to the inverse temperature and the step size in a standard 
Metropolis Monte Carlo programme. If a is too large, the system may stay trapped in a locally favoured 
but globally non-optimal situation. If a is too small, the noise might make the system miss a convenient but 
small structure. The same is true about the noise. Too small a noise makes the programme run very slowly, 
as the clones do not have time to evolve differently, and the selection process is inefficient. Too large a noise 
makes the system miss structures that are dynamically unstable: this is quite crucial when one is dealing 
with highly chaotic systems. 

Last, one of the most difficult parameter to choose is the number of clones used in the simulations [27]. 
Typically, small numbers {N^ ~ 100) will suffice to locate atypical trajectories (say an integrable island in a 
chaotic sea) whereas much larger populations may be needed if one wants to compute accurately the averages 
/i(a,t) or Xi{a,t) (see section 6.2 for a discussion of the convergence of the algorithm as Nc — > oo). 

3. Testing the chaoticity of a system 

An interesting question in dynamical systems concerns the existence of smooth potentials having completely 
chaotic, fully ergodic dynamics. A candidate for this was for many years the Hamiltonian [28] 

n^lipl+pl + x'y') (18) 

The expectation actually turned out not to be justified, as was proven by Dahlqvist & Russberg [29], who 
managed, using an elaborate strategy, to find a tiny (area ~ 10"'') regular island. 

This seems like a good benchmark case to test the Lyapunov Weighted Dynamics [30] . In order to select 
low-Lyapunov trajectories, we use LWD to bias the measure on trajectories by e"'^^*^ with a < and A the 
largest Lyapunov exponent. 

The energy may be taken by rescaling to he E = ^. A small amount of noise is added to the dynamics, 
which rotates the {px,Py) vector by a random angle with typical value ~ V^eSt, where St is the integration 
step, and s controls the amplitude of this energy-conserving noise. As the simulation proceeds, it is slowly 
taken to zero, much in the same way as in standard simulated annealing. Starting from a random point, the 
program falls in a few million steps in the island found by Dahlqvist et Russberg, which shows that it is a 
relatively large target for the programme. We next try to find a new island, but for this we have to avoid 
falling there. This may be easily done by killing all walkers that enter a region including this attractor (for 
example, the rectangle 3.13 < a; < 3.16 and —0.006 < Px ^ 0.004). 

When this is done, the clones find a new island, with area of the order of 10^^". A trajectory within 
this island is depicted in figure 2. One knows that an island has been found by monitoring the finite time 
Lyapunov exponent (cf. figure 3): the beginning of the familiar j regime is a signal that the Lyapunov is 
approaching zero [24]. Note that this strategy works in all dimensionalities, and that it requires no special 
knowledge of the system. Let us note that the noise amplitude (here e = lO^^*') has to be small enough 
that the clones do not "step over" the integrable island in one time-step. The smaller the island one wants 
to detect, the smaller e should be — whence the aforementioned annealing in e done in the simulations. 
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Figure 2. A periodic orbit found by the clones (black curve). The integrable island also contains the image 
of this orbit by the transformation y — —y as well as the rotations of these orbits by 7r/2, tt and 37r/2. The 
trajectory retraces its steps at its ends at the level indicated with dashed lines, where its velocity goes to 
zero. Dashed curves 

orbit found by Dahlqvist &; Russberg [29] . 



^iy^ = 1 enclose the energetically allowed area. The red curve is the shorter periodic 




Figure 3. A logarithmic plot of the finite-time Lyapunov exponent up to time t. The tell-tale sign that a 
region with negligible Lyapunov exponent has been found is the beginning of a 1/4 regime: here it happens 
around t = 1500. We used 1000 clones, cx = -I, e = 10-^°. 



4. Stability of the Lagrange points L4 and L5 



Consider the restricted gravitational three-body problem (see e.g. [31]), with two rescaled masses 1 — /i and 
^ (e.g. Sun and Earth), and a small third body of negligible mass (such as a satellite or asteroid). 

When the large bodies are moving in a circular orbit, there are points where the small mass may be 
placed such that in a rotating frame the three bodies are seen as stationary. These are the so-called Lagrange 
points L1-L5, of which L4 and L5 are stable. When the trajectories of the large bodies are elliptic, with 
an excentricity e > 0, they are not anymore stationary in this rotating frame, but Lagrange points still 
exist. For instance, when /i and e are sufficiently small, it turns out that the vertex of the (time-dependent) 
equilateral triangle defined by the three bodies is a stable trajectory where the small mass moves in phase 
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with the large masses. Many works starting from the one of Danby [32] where devoted to studying this 
situation. 

Denoting the positions of the larger bodies in the rotating frame ri(t) and r2{t), and normalizing lengths 
and times so that the period of a circular orbit is 2tt, the Hamiltonian for the smallest particle in the rotating 
frame is: 

H{r,p)^^{p-Af + V{r) (19) 

with A = {—y,x,0) providing the Coriolis force and V = —^r"^ — (1 — ^)/|r — ri| — — r2\. 

When the orbits of the large bodies are circular, ri{t) = = {—fj,, 0, 0) and r2{t) = rt^ = (1 — /i, 0, 0). 
When e > 0, ri{t) and r2{t) are obtained by solving the Kepler problem; the larger masses librate around 
and rt^ with an amplitude that is larger for larger excentricities. In both cases, the Hamiltonian (19) is 
the only constant of motion. This means that the restricted two-degree of freedom problem, where the small 
particle moves in the same plane as the large ones, is not integrable. It is, in fact, a mixed system with both 
regular and chaotic regions, the former containing the Lagrange points L4 and L5. When the eccentricity e 
increases, or when the masses become comparable, the system is further perturbed by the apparent motion 
in the rotating frame of the large masses. The effect is to reduce further the regular islands. In particular, 
the regular regions containing L4 and L5 shrink and develop resonances, until they eventually disappear. 

We have used the LWD to find the part of the (^, e) plane where regular regions still exist around L4 
and L5. We have used a population of 100 clones and applied a strong bias {a — —10"'), favoring stable 
orbits. We used such a strong bias because as soon as the clones leave the regular vicinity of the Lagrange 
points, they quickly diverge to infinity and have a zero probabilty of coming back. A strong bias thus ensures 
that some clones will always remain within the integrable region. Just as in the previous section, the noise 
required to use LWD is provided by rotating the velocity (p— A) by a weak random angle at each time step, 
thus conserving the Jacobi constant (19). The typical value of this angle is weak: e — 10~^^, and it allows the 
different clones to search for the best situation. For each value of /i, we start all the clones at the Lagrange 
point L4 and adiabatically increase the eccentricity (e(i) = et with e = 10~^ ^ 1), until none of the clones 
is able to remain in the integrable islands and they all diffuse to infinity. This yields the first value Cc at 
which the islands around L4 disappear. We have checked that the results do not depend on the parameter 
e, as long as it is much smaller than all typical time derivatives in the natural dynamics (which are of order 
one) . Note that such a protocol does not give access to the re-entrant part of the stability region (rightmost 
gray part of the figure). In order to obtain the two isolated points on the figure, we slowly increased not 
only e but fj, as well so that these points correspond to the first time the clones leave the integrable island. 
In this way we have obtained automatically (and blindly) the stability region of figure 4. Using a linear 
stability analysis of L4 and L5 , Danby [32] was able to predict a lower bound of this stability region, which 
is very close to the one we observe numerically. This suggests that the regular region surrounding L4 and 
L5 quickly disappears when they become linearly unstable. 



5. The case of integrable dynamical systems 
5.1. Stochastic perturbation of integrable systems 

An integrable system has obviously all its Lyapunov exponents equal to zero. However, a somewhat surprising 
fact is that if we perturb such a system with a weak, additive white noise, it acquires a non-zero Lyapunov 
exponent, defined as the rate of separation of nearby trajectories subjected to the same noise realization. 
Let us consider first a problem with one degree of freedom. For example, for the Hamiltonian 

n = ^+ Viq) (20) 

the equations of motion, once the perturbation ri{t) is added, read: 

q=p 

P=-^+vit) (21) 
aq 
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Figure 4. Stability region (in gray) of the Lagrange points, in terms of eccentricity and mass ratio. 
Boundaries in red are obtained by the LWD while the solid lines are the predictions of Danby [32], based on 
a linear stability analysis of L4 and Lg . The red points connected by a solid line were obtained by increasing 
the eccentricity from e = until none of the clones is able to remain in the integrable island. The two 
disconnected points were obtained by slowly increasing both e and fi so that the first time the clones leave 
the stability region is in the re-entrant region. The number of clones is Nc = 100, the bias is a = -10* and 
the level of noise is set by e = 10"^^. 



Here, ri{t) is a white noise with variance {'r]{t)v{t')) = '2e6{t — t'). For small e, energy drifts slowly, and 
one may prove that there is exponential separation of trajectories along the line of constant energy, induced 
by the noise [33]. The origin of this surprising phenomenon is that noise allows some trajectories to move 
into a faster moving orbit of slightly higher energy, and then come back, a phenomenon that is related to, 
but distinct from, Taylor's diffusion in hydrodynamics [34]. It is, on the other hand, quite different from 
the usual form of 'Noise Induced Chaos', which is originated by the noise-induced occasional visits a system 
makes of unstable chaotic regions [35], there being in our case no chaotic regions at all. 

Clearly, for a harmonic oscillator the present mechanism cannot work, as all orbits have the same period, 
regardless of the energy. As we will show elsewhere [33], when the variance of the noise e — 0, the diffusion 
occurs mostly tangentially to the tori on time-scales up to t e~^^^. The induced divergence of two nearby 
trajectories is exponential and the corresponding Lyapunov exponent can be computed explicitly [33]: 

where / and (f> are the action and angle variables of the problem, u = (j) = ^ is the frequency. The variable 
q{t) is, in the absence of noise, a periodic function of 4>, the overbar in (22) denotes average over a cycle of 
the noiseless dynamics. The Lyapunov time is proportional to e~3^ while the diffusion away from a torus 
takes times of order e~^: for small e there is then a wide range of times ^ t <C where the definition 
of a Lyapunov exponent associated with a torus makes sense. 

An interesting limit of equation (22) corresponds to the approach of a separatrix of energy E^, where 
the frequency uj vanishes. One may easily check that, quite generally, in terms of AE — E ~ E^: 

u.iAE^O)^-^^0 (23) 
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so that: 



■ cx 



1 



dp "~ dE AE\ log AE\3 ^ °° ^^"^^ 

Hence, as one may have expected, the Lyapunov exponent becomes large as one approaches a separatrix, on 
which it becomes larger than ©(e^/^). 
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Figure 5. LWD with 2 000 clones for a system defined by (25) and e = 10~®. Top: «i = 1 and cji>2 = 0. 
Position of the clones for t ~ 5 000 (dark green), t ^ 10 000 (purple) and t ~ 15 000 (gray). Bottom: 
ai 2 = 1 and 03,4 = 0. Stationary positions of the clones (here t ~ 2 000 in gray). 



5.2. The special case a = 1 

In order to see what kind of structures dominates the dynamical partition function Zp = (^e^^i+---+^p)t'j ^ 
let us consider a system of two degrees of freedom that consists of a cartesian product of two systems of 
one degree of freedom. In figure 5 we show the result of LWD for a cartesian product of two double wells 
T~1-{q^p) — 2Pi/'^ + (if ~ 1)^/4' with additive Gaussian white noise of variance 2e: 

Qi ^Pi Pi = -qiiqf -l) + V2eT], (25) 
This system has two normally hyperbolic invariant manifold (NHIM) with one unstable direction defined by 
Qi = Pi = for the first one and 92 = P2 = for the second one. They correspond to the cartesian products 
between the flat measure § over one double well and the saddle point of the other double well. It also has 
one NHIM with two unstable directions defined hy qi — q2 — pi = P2 = which is the cartesian product of 
the two saddle-points. The walkers are biased to search for trajectories with either an atypically large Ai 
(tti — 1 and ai>2 = 0) or an atypically large sum Ai -I- A2 (ai = a2 — 1, ct^ = a4 = 0). 

§ Note that (25) is the limit of a Kramers equation with friction jp and temperature kT when 7 — and kT —> 00 with 
e = ■ykT constant. Without cloning, the steady-state measure is thus the infinite temperature limit of a Boltzmann weight, i.e. 
the flat measure. 
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When we bias the walkers with e'*'^*, we stabilize the unstable manifold of one NHIM with one unstable 
direction: the clones are localized on the separatrix in one of the double wells while they diffuse freely in 
the other one (since the noise does not conserve the total energy). On the other hand, e^'*''^"'"'^^^* stabilizes 
the unstable manifold of the NHIM with two unstable directions, i.e. the cartesian product of the two 
separatrices. 

5.3. The most regular trajectories (a <Q) 

As already mentioned above, choosing a = 1 is such that the walkers populate the separatrix uniformly: the 
stretching and compressing in times of acceleration and deceleration are exactly compensated by the cloning 
rate [21]. In the following we consider the converse case with a < 0. 

In a Hamiltonian system, the motion in a regular region is constrained to a torus, labeled by the 
constants of motion la, and spanned by the dynamic curves 9a = §y^. By changing the values of the la, 
one moves to different nested tori. Consider for definiteness a system with two degrees of freedom with two 
constants of motion Ji, I2 and the corresponding angles 61, 62 which span the two-dimensional surface of the 
torus. The innermost of these tori reduces to a ring, spanned by only one of the two angles, say 9i. In the 
vicinity of this ring, the dynamics of /2,^2 generically reduces to a harmonic oscillator. In mixed systems, 
such a degenerate trajectory corresponds to the center of an integrable island. 

Let us now show how the LWD can be used to find such a degenerate structure. As shown in section 5.1, 
the noisc-induccd Lyapunov exponent in integrable systems come from the anharmonicities that make 
7i"{I) ^ 0. This may be seen in equation (22) by checking that in the case of a harmonic oscillator, 
for which % — I , the Lyapunov exponent vanishes. Looking for the trajectory with the smallest Lyapunov 
exponent in the presence of noise, amounts in fact to looking for the "most harmonic" of trajectories. In the 
case of nested tori, this generically corresponds to the innermost trajectory, which is degenerate. 

As an example, consider the restricted three-body problem with /i = 0.1 and the large masses performing 
circular orbits. As mentioned above, the system has many regular islands, that correspond to non-chaotic 
orbits, many of which have been classified [36]. Figure 6 concerns one such regular island: the typical 
trajectories in there, such as the red one, are rather complicated in real space. Now consider the degenerate 
"innermost" torus in the island, in black: it corresponds to the "flower" structure which is stationary in 
the rotating frame. Indeed, we recognize every other regular trajectories of the island as essentially this 
trajectory, with, in addition, librations of this basic pattern. Clearly, even for analytic purposes, it seems 
useful to be able to locate such a structure. (Their role is, by the way, reminiscent of the role played 
by "inherent structures" in glassy problems where all the vibrational motion is eliminated by treating as 
single configuration all points belonging to a same basin of attraction of the zero-temperature dynamics). 
Furthermore, while it is possible to "guess" that the red trajectory of figure 6 corresponds to librations 
around the black one, this would be impossible in higher dimensions while the LWD would nevertheless 
have no trouble finding the proper degenerate trajectory. Figure 7 shows, for the example above, the clones 
migrating to the center of an integrable island, when selected for low Lyapunov exponent. 

Higher dimensional structures may be selected for an integrable problem with more degrees of freedom. 
For instance, minimizing the sum of the smallest p positive exponents Ai -I- ... -t- Ap leads to tori that are p 
times degenerate (i.e. they are N — p dimensional) . These considerations might be interesting to apply to 
extended integrable systems, such as the Toda lattice, because it might give a method to construct soliton 
solutions. 
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Figure 6. Simulation of the restricted three-body problem (19) with an eccentricity e = and fi = 0.1. 
Red: quasi-periodic stable orbit. Black; the periodic orbit at the center of the corresponding island of 
stability. 
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Figure 7. Poincare section {y = 0, E = 'H(L2)) corresponding to figure 6. The clones are drifting from 
the initial position (a) to a nearby integrablc island (b corresponds to the red trajectory of fig 6). They 
then migrate to the center of the island (c corresponds to the black trajectory of figure 6). Here H is given 
by (19), we used 100 clones, a = — lO", and a noise intensity e = 10~^. The color code shows the time spent 
since the beginning of the simulation. 
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6. Spatially extended systems 



6.1. Fluctuations of several Lyapunov exponents: the Fermi- Pasta-Ulam chain 

Let us now turn towards spatially extended systems and enter the domain of condensed matter by considering 
the Fermi-Pasta-Ulam-Tsingou chain defined by the Hamiltonian 



with xi^^i — xi. This corresponds to a chain of L particles coupled with anharmonic springs. In the p — Q 
limit, the system is integrable and the chaoticity increases as /3 and % increase [37]. In this paper, we consider 
the case /3 = 0.1 and H = L. In previous studies [24, 25] we showed that biasing this system in favor of 
chaotic trajectories reveals chaotic breathers, whereas regular trajectories correspond to a gas of solitons 
propagating ballistically through the system. The unbiased case corresponds to a mixture of short-lived 
solitons and breathers superimposed with thermal fluctuations. 

Two questions arise naturally when considering spatially extended systems. First, what happens when 
we consider Lyapunov exponents beyond the largest? Second, how do the fluctuations scale with the system 
size? Answering the latter in the FPU case is a difficult task since even the average values already converge 
very slowly with the system size [38]. Below, we thus address the former question. When using the LWD 
to sample the large deviations of the Hamiltonian dynamics associated with (26), we use the energy- and 
impulsion-conserving noise described in Appendix A. Since the chaoticity of the FPU chain depends on the 
energy of the oscillators, it is indeed quite important to prevent the system from decreasing its chaoticity by 
either decreasing its total energy or by using all its energy to create a uniform rotation of all the oscillators. 

In [7] it was shown that above a certain energy, a modulational instability of the FPU chain leads to 
the formation of a chaotic breather: the energy condensates on a small number of degrees of freedom and 
the largest Lyapunov exponent increases substantially. In [24] we showed that these chaotic breathers are 
the trajectories that contribute the most to Z{a,t) for large positive a. 

As we can see in figure 8, the appearance of a chaotic breather creates a gap in the Lyapunov spectrum, 
with the largest Lyapunov exponent increasing by a factor 5, and the second one by a factor 2. The plot of 
figure 8 includes cloning and noise, so that one may think that the individual trajectories differ completely 
from the underlying FPU dynamics, but this is not so: using the positions and momenta of a given clone 
as an initial condition for a numerical simulation without cloning and noise shows that the breathers are 
indeed a proper solution of the FPU chain. As explained in the introduction, the cloning simply stabilizes a 
metastable solution of the equations of motion. Since clones are constantly killed or copied as time goes one, 
one could be surprised that the trajectory in figure 8 does not show clear cloning events where, for instance, 
the breather would suddenly disappear and reappear few hundreds of site further. This can happen but is 
quite rare because, once a breather has been found, all the clones in a small population tend to be very 
similar. This "degeneracy" of the clone population, that sample the vicinity of one interesting trajectory 
while many exist, is one of the reason why in section 6.2 we will need much larger population to accurately 
compute the average Z{a,t) = (6""^^*). 

The few first tangent vectors are localized on the breather, and one can thus wonder what happens when 
we require that other Lyapunov exponents be large. Either some internal degrees of freedom of the breather 
can be excited, leading to an increase in Lyapunov exponents beyond Ai, or the system may divide the 
energy of the breather into several breathers, each of them corresponding to one atypically large Lyapunov 
exponent. A trajectory containing k breathers is thus a natural candidate for large deviations of the k first 
Lyapunov exponents but, as was shown in [7], the breathers have a tendency to merge upon collision which 
makes multi-breathers solution unstable. 

As we show here, the breather-breather interaction is actually more complicated and multi-breather 
solutions indeed exist. In figure 9 we show that biasing with ai — ai>3 = and ck2 > favors the 
appearance of two breathers. Another Lyapunov exponent, A2, increases and the gap is now between A2 and 
A3. The same is true for A3 and A4 that generate solutions with 3 or 4 breathers (see figure 9). A large A2 
also implies a large Ai since Ai > A2 by construction. Using ai > and a2 > can thus also produce two 
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Figure 8. Simulation of a FPU chain with the Lyapunov Weighted Dynamics and parameters L = 512, 
a = 5L, 200 clones, dt = 10"'^, e = 10"'^. Left: Energy of each site (for one clone) as a function of time. 
Centre: vs i for i S [1, 10] at t = 10* for the biased (blue) and unbiased (red) case. The largest and 
second largest exponents are clearly larger in the presence of a breather. The inset shows Xi{t) for t < 
for i 6 [1, 10]. Right: v^^ + v^^ as a function of time and lattice site i for the first tangent vector, which is 
shown to be localized along the breather. 
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Figure 9. Multi-breather solutions in a L = 512 system obtained using LWD with 200 clones and parameters 
dt = 10~^ and e = 10~^. Left: ai^2 = 0, 02 = 5L yields two breathers. Center: Qi^^a = 0, 03 = 5L 
yields three breathers. Right: a^^^ = 0, 04 = 5L yields four breathers. 



breathers. However, as can be seen in figure 8, a single breather aheady has a larger A2 so that either one- 
or two-breather solutions dominate the biased measure depending on the value of ai and a2- In practice, to 
maintain a small cloning rate and nevertheless be sure of seeing two breathers, we used ai = and increased 

Q!2. 

The question as to how one transitions from typical profiles (homogeneous for ai ~ 0) to atypical ones 
(with breathers for large enough ai) is an interesting one. In particular, it would be interesting to know 
whether there is a phase transition for critical values and if so how the a^s scale with the system size. 
This is left for future work and we shall only here prove that the method can indeed single out atypical 
trajectories in spatially-extended systems. 

Interestingly, it was reported in [7] that the formation of a single breather starting from a short 
wavelength perturbation of the homogeneous state relied on the merging of breathers upon collisions. The 
only way to observe a multi-breather solution thus seemed to have a system large enough that two breathers 
would not have the time to meet during their life time. If that was the case, the solution presented in 
figure 9 would thus be artefact due to noise and cloning. To check that this is not the case, we pursued 
the simulations without noise and cloning, and found that the breathers indeed remain (meta)stable, even 
after encounters. This revealed a great diversity in breather-breather interactions, as plotted in figure 10. 
In addition to the coalescence reported in [7], the breather-breather interaction can also be repulsive, the 
breathers turning back before colliding, or they can simply cross each other with little energy exchange. 

To conclude, multi-breather solutions thus exist and were not detected previously in the literature due 
to the difficulty of finding the proper initial conditions. This is here unessential thanks to the clone dynamics 
which automatically locate these highly aptypical trajectories ||. 

II We thank the anonymous referee for suggesting this synthetic summary. 
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Figure 10. Close-up on two types of breather-breather interactions: crossing and repulsion, from left to 
right. We used the position and impulsion of a given clone as initial conditions for standard Hamiltonian 
dynamics. 



6.2. Large deviation functions: the example of coupled maps 

Beyond the detection of rare trajectories, one may be interested in measuring their actual probabihty of 
occurrence. Many systems exhibit intermittency, and show an apparent metastabihty in space-time of 
trajectories of different chaoticity levels. The construction of the cumulant generating function gives access 
to a free-energy like quantity that allows one to export the language of phase transitions to such dynamical 
phase coexistence [10]. For instance, the average Lyapunov exponent Ai(a) in the clone population plays 
the role of an order parameter since it is the first derivative of the dynamical free energy /x(q;): 

(27) 

As far as we know, the cumulant generating function of Lyapunov exponents has only be computed 
for low dimensional models, such as maps of the interval [39], with the notable exception of the Lorcntz 
gas [40, 41, 42]. Recently, Vanneste extended the algorithms introduced in [24] to compute the generating 
function of the Lyapunov exponents of products of random matrices [43] . He showed how the cloning method 
helps reducing the variance of estimators based on finite number of clones and computed the large deviation 
functions for products of small random matrices (up to 8 x 8). Here we show how the LWD can be used 
to compute the cumulant generating function of a spatially extended dynamical system: a one-dimensional 
lattice of coupled skewed tent maps (CSTM) [44]. 

We consider L real variables Xi e [0, 1] whose discrete-time dynamics are given by 

bx if a; < T 

x,(i+l) =/(x,(i))+^[/(x,+i(i))+/(x,_i(i))-2/(x,(i))] f{x)^{^_^ (28) 

otherwise 



1 



1 



b 

When D = 0, all the maps are uncoupled and fully chaotic. The stationary measure is uniform and all 
Lyapunov exponents are equal to [17] 



A = log \f{x)\dx = ^ logfe + (^1 - log ^ 



(29) 



When D is increased, the maps become locally coupled and the dynamics are made more regular by the 
diffusion terms, whence a decrease of the Lyapunov exponents (we use D = 0.1 and 6 = 4 in the following). 
This system is a good test for our method because it is simpler to simulate than continuous time dynamics 
and its Lyapunov spectrum converges quite rapidly when L increases (see figure 11), compared to more 
complicated systems such as the FPU chain of section 6.1. Furthermore, we may check whether our algorithm 
allows us to go beyond the existing works on the fiuctuations of Lyapunov exponents in coupled-map lattices, 
which have been limited to the study of the Gaussian regime [15]. 
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Figure 11. Spectrum of Lyapunov exponents for L = 20, 40, 80, 160, 320 and 640. The exponents are 
averaged over a large number of runs, from 5 000 for L = 640 to 500 000 for L = 20. For practical purpose, 
the spectrum has converged when L > 40. 



To compute the Lyapunov exponents, we need to linearize the dynamics (28) to determine the time 
evolution of a tangent vectors u: 

u,{t + 1) = (1 - 2D)f{x,{t))u,{t) + D[f\x,+{)u,+i{t) + (30) 
b if 2; < ^ 

/'(^) = <! ^ ^ (31) 

otherwise 



1 - h 

Again, to use the LWD we will need to make (28) stochastic. To do so, we use the shift 

elJ 

Xi{t) = x^{t) + ^ X min(xi(i), 1 - x^{t)) (32) 

where J7 is a random number sampled uniformly in [—1, 1], e set the scale of the noise and min(xi, 1 — xi) 
makes sure Xi{t) remains in [0,1] for e < 2. We then replace Xi{t) by Xi{t) in (28). The "noise" (32) is 
clearly not the Gaussian white noise discussed in section 7 for Hamiltonian dynamics. One thus has to be 
careful and check that the observed fluctuations of Lyapunov exponents do not depend on the details and 
intensity of the noise (see below). 

We first use brute-force simulations for L = 40 to estimate the Gaussian scaling of the fluctuations of 
Ai and check the effect of the noise (32) on the fluctuations of the coupled maps. After an initial run of 
500 time-steps, the Lyapunov exponents are computed over the next t = 10^ iterations of the maps. For 
N = 10® simulations, figure 12 shows that the pdf is well approximated by a Gaussian P(3(Ai,t): 



P,i„,(Ai,i) ~ PG(Ai,t) = with (Ai)~ 0.372 ~ 2.03 x 10-^(33) 



Note that the standard deviation of Ai is o- = (j/vi ^ 4. 10^ (Ai), showing that these simulations do 
not allow to sample very far from the average value. Nevertheless, figure 12 shows the beginning of the large 
deviations regime since there seems to be a small but systematic asymmetric deviations from the Gaussian 
approximation when P{\i,t) ~ 10"** — 10^®. Importantly, simulations with levels of noise e ~ 10^^ and 
e = 10^^ show no difference in P(Ai,t), so that e = 10^^ will thus be used in the following. 
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Figure 12. Left: Psim.(^l)t) constructed from direct sampling of 10® simulations for L = 40, D = 0.1, 
fe = 4, £ = 10~^ (red) or £ = 10"* (blue). The black line is the Gaussian approximation. Right: same plot 
in logscale. One sees a small asymmetric deviation from the Gaussian when P(Ai,f) ^ 10^* — 10^®. 



From the Gaussian approximation, one can construct a quadratic approximation to n{a) 



1 



t 



1 



Ai^(a) = -log / e"^i*FG(Ai,t)dAi =a{\i) + -a'a' (34) 



2 



Similarly, the value of the Lyapunov exponent \i{a) that dominates the biased measure for a fixed value of 
a is approximated by 



(Aie"^i* 

/gCiAit 



^1 (") = / = t) = (Ai) + aa' (35) 



One can then compare the Gaussian approximations with /^(a) and Ai(a) obtained using the LWD of 
section 7.3. The simulations in figure 13 were obtained using varying numbers of clones, from Nc — 100 to 
= 160 000. For each value of a and N^, Ai(a) and /i(a) were obtained by averaging over 10 runs of N^. 
clones. For a < the convergence with the number of clones is very fast, the difference in Ai(a) between 
Nc = 100 and Nc = 160 000 for a = —3 is less than 0.5%. For a > 0, the convergence with Nc is slower 
when a gets larger. To estimate the speed of convergence, one can check for which value of a the difference 
between the estimators of Ai(a) for Nc and 160 000 clones goes above 5%. One gets a = 1.4, 2.0, 2.2, 2.4 
for Nc = 100, 500, 1 000 and iV^ = 2 000. The estimators for N^ = 10 000, 20 000, 40 000 and 80 000 always 
remain within 3.5%, 2%, 1% and 0.5% of the values obtained for N^ = 160 000. The same simulations done 
with different noise intensities (e = 10^'' or e = 5. 10^^) lead to the same asymptotic values for Ai(a) and 
/i(a) when Nc — )• oo, though the convergence rate is slower at smaller noise intensities (see figure 14). For 
Nc = 160 000, the difference between Ai(a) for the three noise intensities is smaller than 0.5% for a < 2.2 
whereas it goes up to 2% at a = 3.0. 

Let us highlight here that the number of clones used to compute Ai(a) is much larger than what was 
used in section 6.1 to detect the breathers. This is mostly due to the fact that in section 6.1, we simply 
wanted to detect an atypical trajectory, without quantifying its rarity, whereas in this section we aim at 
computing averages that require sampling over many different trajectories. There is no obvious rule as to 
when a large number of clones will be needed to compute Ai(a) (for a < the convergence is very fast, 
for instance). Also, many tricks exist in the literature for Sequential Monte Carlo simulations (see [45] and 
reference therein) that could/should probably be extended to the LWD to speed up the convergence of the 
algorithm. 

In figure 13 we see that the LWD indeed samples far beyond the Gaussian regime: the values of Xi{a) 
vary over an interval of order (Ai) whereas a was less than 1% of (Ai). Let us also note that the deviations 
from the Gaussian scaling are more clearly seen in Ai(q;) than in /Lt(a). This can be understood because 
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Figure 14. Convergence of Ai(q) as a function of the clone number Nc for three levels of noise. A given 
color corresponds to a given value of a virhile different symbols correspond to different levels of noise. 



the latter is nothing but the integral of the former. Since A(a) is quite small (~ 0.5 here), we need a large 
interval of a for the error on Ai(a) to have an impact on fi(a). Nevertheless, it is quite reassuring that our 
simulations agree with the Gaussian scaling when a ~ 0. 

Since there are neither analytical nor numerical benchmarks for our computations of fi{a) and Ai(a), 
we compared the results obtained using two different types of cloning algorithms. We used both the 
methods described in the core and at the end of section 7.3, that use resampling of the clone population 
respectively clone-by-clone or globally. Both methods gave similar results, the clone-by-clone method (which 
we generically refer to when we speak about LWD) being substantially faster that the global resampling 
scheme. We have thus shown that LWD can be used in spatially extended systems to compute the cumulant 
generating functions beyond the Gaussian regime. This now opens the way to systematic studies of the 
dependence of /x(a) and Ai(a) with the system size, which should allow us to discuss the possibihty of 
dynamical phase transitions separating different scalings with the system size of the exponents fluctuations. 
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7. Method 



7.1. Definition of finite-time Lyapunov exponents 

Let us consider a system evolving under Hamiltonian dynamics, which can be written in a concise form 

{x = {qi,...,qN,Pi,---,PN) 
^ an &}£ _dn (36) 
dpi''"' dpN ' dqi''"' dqN 

The Lyapunov exponents characterize the exponential convergence or divergence between two nearby 
trajectories. They are usually defined by considering the linearized dynamics of a tangent vectors uj: 

u:^-Au: with A„- - -^4!^ (37) 

dx, ^ ' 

There are two ways of defining finite-time Lyapunov exponents. For a given trajectory, let us call U{t) 
the matrix solution of the linear equation U{t) = — A[x(t)]C/(t) and Ai{t) > A2{t) . . . > A^it) the N first 
eigenvalues of U^U. We then obtain a first definition of the finite-time Lyapunov exponents as 

A.(t)- ^logA(i) (38) 
Alternatively, one can look at a /c- volume Vk{t) = uji{t) A u)2{t) . . . A oJk{t) and define 

\k{t) = lim -\og\Vk{t)\ ~ lim ^log|14-i(i)l (39) 

which thus corresponds to the exponential growth in the fcth dimension of the fc- volume 14 . In the limit 
t — oo, for generic choices of the vectors both definitions coincide, although they generically differ at 
finite time. Expression (38) is more satisfactory because the exponents do not depend on the choice of 
the ijJiS, but (39) is much simpler to implement numerically. The method developed in this paper uses a 
Gram-Schmidt orthogonalization procedure that relies on definition (39), but it would be interesting to see 
how the fluctuations of finite-time Lyapunov exponents differ using definition (38). A natural way to do so 
may be to use the recent algorithms introduced to compute covariant Lyapunov vectors [44, 46]. 

Let us note furthermore that both (38) and (39) explicitly depend on the choice of trajectory x{t) around 
which the linearized dynamics (37) is studied, and hence we could/should write Ai[t, x(0)]. For non-ergodic 
systems, even the long time limit \i{t — >■ oo) still depends on x(0) whereas for ergodic systems Osedelets 
ergodic multiplicative theorem ensures that \i{t — >■ oo) adopts the same value for (almost all) xq [18]. In 
the following we drop the dependence on xq but one should remember that it is an important source of 
fluctuations of Xi{t). 

Let us show how to relate equation (39) to the tangent dynamics (37). To do so, starting from the 
vector ijJi, we obtain a set of orthogonal vectors Ui which generate the same oriented volume Vk through 



i=i ^ ' 
whose dynamics can be inferred from (37) 

4-^ Ui ■ Aui + Ui ■ Am," 
w, = -A«, + 5]«,^ f-^ ^ (41) 



\Ui, 



One directly checks that ^^(Mi • Mj) = so that the UiS remain orthogonal if they were chosen so at t = 0. 
Furthermore, A u;2 • • ■ A — Ui A U2 ■ . ■ A u^., since an exterior product involving twice the same vector 
vanishes. 

Let us then introduce the inner product 

{ui A U2 . ■ ■ A Uk\vi A V2 ■ . . A Vk) = det M with Mij = Ui ■ Vj (42) 



19 



to express the norm of Vk through |Vfep = {ui{t) Au2{t) . . . Auk{t)\ui{t) Au2{t) . . . A Uk{t)) so that the thue 
evolution of Vk is given by 

^iVfeP = 2^(mi a . . . a ttj a . . . a Mfeliti a . . . Mfe) (43) 

j 

All terms in iij parallel to Ui^j vanish, yielding 

^|VfcP = -2^(mi a . .. a Amj a ... AMfeltti A...Uk) 

3 

We thus have to compute a sum of determinants of the type 



(44) 



l^^il 






\U2\ 







Ui ■ Auj U2 ■ Auj 




Uj ■ Auj 






Uk ■ Auk 


\Uk\^ 



= l«i| 



• AMj|Mj + iP ...\Uk\ 



and thus 



dt 



whose evolution read 



1^41' = -2 Jl • Av, ^ -2|t4p ^t;, • Av, 

i—l i—1 i—1 

where we have introduced the unitary vectors Vi 

Vi = -Avi + Vi{vi ■ Avi) + ^ Vj{vi ■ Avj + vj ■ Avi) 

Finally (45) can be solved to yield 

\Vkit)\=e-^o^'{^L^-^■^-^}\VkiO)\ 
One can thus rewrite (39) as 



1 /•* 

- / dt{vk- Avk} 
^ Jo 



(45) 



(46) 



(47) 



(48) 



This formalism may seem complicated but (4G) simply yields the continuous time evolution of the Gram- 
Schmidt vectors [47], while the Afc(t) defined in (48) can be computed using the rescaling factors that appear 
in the renormalization steps of the Gram-Schmidt procedure. This procedure indeed corresponds to evolving 
k unitary vectors Vi initially chosen orthogonal with the tangent dynamics 

v^ = -Av, (49) 

At every time step t — rdt we re-orthogonalize them and normalize them, calling Sj (rdt) the rescaling factor 
of Vj. Then, when dt — ^ 0, the dynamics of Vi tends toward (46) and 

t/dt k 

' ' " ] V — ^ ^^^^ 



t/dt k k k 

n n *J exp ( - / dt ^ Vj ■ Avj^ ^ exp {t ^ (t) 



Let us now show how this formalism allows us to sample the large deviations of Lyapunov exponents. 
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1.2. Large deviation formalism 



In general, one needs different levels of control of the details of the algorithm depending on whether one 
wants to compute a dynamical free energy (or topological pressure) or simply look for atypical trajectories. 
One important point is that for the simulations, it is always necessary to add a small amount of noise to the 
system. We shall assume that this is the case below, but one should remember that in some cases the noise 
is simply an artifact and has to be chosen as small as possible. We thus consider a slightly different version 
of the dynamics (36): 

9'H ■ dn ,— 
dp, dq, 

where jy, are unitary Gaussian white noises. To compute the dynamical partition function 

Z(a,t) = (e"^^*) (52) 

we first need to define the average (...)• In deterministic dynamical systems, finding a putative stationary 
measure is difficult and the addition of a Gaussian white noise simplifies drastically the situation. For 
instance, (51) would lead to a flat measure over the whole phase space, while considering a stochastic force 
acting tangentially to the energy surface would make the microcanonical measure be the correct steady- 
state (see Appendix A). In the latter case, the average (52) then amounts to choosing the initial conditions 
uniformly over the energy surface and averaging over the noise realizations. 

From the mathematical point of view, the addition of a small noise changes the nature of the system. 
The question as to whether one recovers, in the small noise limit, the steady-state measure of the underlying 
deterministic system dates back to Kolmogorov (see [48]). This question of "stochastic stability" [49] is 
actually a natural way for a physicist to define the steady-state measure of a dynamical system, and it is 
thus fortunate that when SRB measures exist, they can indeed be recovered as small noise limits of stochastic 
dynamical systems (see [49] for a proper mathematical presentation of this procedure). 



7.2.1. Largest Lyapunov exponent For simplicity, we first consider the large deviations of the largest 
Lyapunov exponent Ai and later generalize the formalism to fluctuations of several exponents. Furthermore, 
we will use the Gaussian noises introduced in equation (51) and discuss in Appendix A how to proceed with 
conserved noises. The probability density obeys a modified Liouville equation \. 

^ -H^P[q,P,t) He=^E-^ + ^^-^^ (53) 

Including the dynamics of the tangent vector v, and using the notation x = {q,p), this becomes 
dPix,v,t) ^ d 



dt 



= -HP{x, v,t) = -(^H,-Y,g^[Yl ^'J^J + ) Pi^' t) (54) 



i=l j 



where we have introduced N{v) ~ ^Y^'ijLi''^i^ij^j fo^' clarity. We can then define a distribution of Xi{t) 
through 

P{Xi,t) ^ J 'D[x,v]dxodvo5[Xi - Xi{x,v,t)]P{xo,v^) (55) 

where Xi{x,v,t) is defined in (48). The evolution of the unitary tangent vectors v is given in (46). 

We will now rewrite (52) as the path integral of a generalized evolution operator in which we will then 
read the algorithm used to bias the trajectories with a weight exp(aAii). We wish to compute 

(e"^i*) = (e-"^'J- ^0 ^^^'j^jdt^ (56) 

^ Derivatives apply to everything on their right when they appear as operator, as in and not when explicitly applied to a 
function, as in S^. 

' ox 
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which can be rewritten as 

2N 2 

«=1 j kl 

We thus integrate over all possible trajectories starting from all possible initial conditions. The 5 functions 
constrain these trajectories to be solutions of the equations of motion and the Gaussian weight is the 
probability of a given noise realization. Since we will only use this path integral to recognize an evolution 
operator, we do not need to worry too much about operator ordering and time-discretization of multiplicative 
noises. Using the complex representation of 5 functions and doing explicitly the Gaussian integral over the 
noise, we get: 

V[x, V, X, v] da;odt;o exp J dt{ E [mi + PiPi + ePi ~ + 'dq-^^^ ^^^^ 

2N 2N 

+ ^[viVi+v,{'^kijVj+ViN{v))]+aN{v)] P{xo,Vo) (59) 

where x, v are imaginary fields. This path integral corresponds to the matrix element 

(e«At^ ^ (-|e-*(^-«^(''))|Po(xo,t;o)) (60) 
where H is given by 

2^ , 2N 

^ = - E [ E ^ij^'j + ^(^)^*] (61) 

The cumulant generating function /x(a) = | logZ(Q!,t) is then given by minus the smallest eigenvalue 
oi H — aN{v). The evolution operator H — aN(v) docs not correspond to a standard Langevin equation 
since it does not conserve the total probability. The evolution equation 

P = -{H - aN{v))P (62) 

indeed implies that 

y dxdvP{x,v;t) = —a J dxdv '^^ViAijVjP{x,v;t) (63) 

which does not generically vanish. This means that if one wishes to represent P(x,v,t) by a population of 
points in the space {x, v), the number of points is not conserved. Each copy the system is thus replicated (or 
pruned) at a rate aN{v), i.e. a times the stretching rate of the tangent vector u{t). Before we turn to the 
numerical implementation of this dynamics, let us show how the discussion above extends to the fluctuations 
of several Lyapunov exponents. 

7.2.2. Large deviations of several Lyapunov exponents Let us consider the joint fluctuations of the L flrst 
Lyapunov exponents. We thus wish to compute the generating function 

Z{cx, t) = (e- ^^=1 \ (64) 



where a = {ai. . . . ,aL) is a vector of "biases" corresponding to each Lyapunov exponents. Following a 
similar path as in the previous subsection, we can write Z as 

Z{oc, t) = (-|e-*(^-^^- |Po(9o,Po, ^^1,0, • • • , vl,o)) ~ e*''(°) (65) 
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where //(a) is now minus the smaUest eigenvalue of the evolution operator H — X]fc=i oikN{vk), where 



N 



^ 1 dp; dpi dqi dqi dpi I 



1=1 

L 



A^K) - -t^fc • A-yfc (68) 

We thus have to represent P{x, Vi, . . . , v^; t) by a population of walkers in the space {x, Vi, . . . , v^). The 
pruning or cloning rate is now given by the product of each of the factors akN{vk), which are nothing but 
ak times the stretching rate of the vector Vk- 

7.3. Algorithm: the Lyapunov Weighted Dynamics (LWD) 

Let us now turn to the numerical implementation of the dynamics generated by J? — J2k'^kN{vk). To 
do so, we will use a population dynamics that allows one to simulate dynamics that do not conserve 
probability. This type of algorithm is closely related to Sequential Monte Carlo dynamics or Diffusion 
Monte Carlo [50, 51, 52, 53, 54, 55, 56]. An alternative to this type of approach is to carry Metropolis like 
Monte Carlo directly in the trajectory space, as does Transition path sampling [57] which has been used to 
compute large deviation functions of dynamical observables in models of glass forming liquids [58] (see [59] 
for a version of transition path sampling that uses Lyapunov exponents to weigh the trajectories). 

We consider a population of Nc copies of the system with positions and momenta q and p. To each 
copy j, we attach L tangent vectors Vi. We then choose a time-step dt and evolve the system over T time 
steps so that the total time is t = Tdt. For t ~ 0, the Nc copies start from arbitrary initial conditions. At 
each time step t' = ndt: 

\Y] For each copy j 

• {QtP) evolves with the noisy Hamiltonian dynamics (51), 

• each tangent vector Vi evolves with the tangent dynamics (49) 

• we use a Gram-Schmidt procedure to turn the VjS back into an orthonormal basis: for each i, 

* we make Vi orthogonal to each Vk<,i'. Vi ^ Vi — X]fe=i '^ki'^i ■ v^) 

* we define Si(n) = \vi\ and then normalize Vi: Vi jrjy^Vi 

• we use the rescaling factors of the tangent vectors to compute the weight Wj{n) = Si{n)°'' 

|~2^ we then compute the average weight 

^(") = ]^E^jW (69) 

[3] each copy is then replaced on average by Wj{n)/ R{n) copies. To do so, we pull a random number Sj 
between and 1. The copy j is replaced by+ r — [ej + Wj{n)/R{n)\, 

• if T = 0, we delete the copy 

• if T > 1, we create r — 1 new identical copies, 

pl] at this stage, we have — J2j l^j + Wj{n)/R{n)\ copies. On average, we are left with Nc clones since 
(Nc) — Nc- However, we do not want the number of clones to diffuse and thus keep it strictly constant: 
we kill at random Nc — Nc clones if Nc > Nc and we duplicate at random Nc — Nc otherwise. 

The dynamical partition function and free energy are then given by 

T T 

Z{a,t)^l[R{n) ^i{a,t)^-J2^ogR{n) (70) 

n— 1 n—1 
[xj is the largest integer smaller than x. 
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The distribution of trajectories described by the walkers converge to the natural distribution of the dynamical 
systems, weighted with an extra factor exp(^j, afeAfct). Note that as for any Monte Carlo simulations, 
metastability can be a problem. The clone population can indeed get trapped in a metastable "state" in the 
trajectory space composed of trajectories that are locally favored but do not globally dominate the average 
Z{a) = (6^'="'=^'=*). In such cases, it is better to run several simulations in parallel and average /i(a, t) over 
these simulations rather than increase the number of clones in one simulation. 

Global resampling scheme. An alternative to step [2] — [4^ is to resample the whole population according to 
the weights WjS. This can be done using tower sampling [60]: we construct the cumulative weight function 
C(0) — wq and C{j > 1) = C{j — 1) + Wj. Then we pull Nc random numbers between and C{Nc). Every 
time such a number falls in the interval [C{j), C{j + 1)] we make a new copy of the clone j. This strategy 
is often used in the Sequential Monte Carlo literature [45]. We tried both strategies which gave similar 
results. For a ~ 0, the LWD presented above barely requires any cloning since r ~ 1 for all clones. On the 
other hand, the global resampling method amounts for a ~ to choosing Nc random number between and 
Nc- On average, there is one number in each interval [n,n + 1] but because of fluctuations many intervals 
will contain more than one random number and many intervals will contain no random number. The global 
resampling then typically involve many cloning events (typically N^/e for a = and N,, 3> 1 *). Fluctuations 
are thus stronger for the global resampling strategy and the overhead due to the cloning is typically larger, 
whence a slower algorithm. The code is, however, much simpler to write and a large literature exists on how 
to implement it efficiently (see [45] and reference therein). 

Another alternative to the LWD, which constructs a canonical measure over the trajectory space using 
an inverse temperature a, would be to try and implement a multi-canonical algorithm — a strategy that was 
for instance followed in [61]. Again, this is left for future work. 

8. Conclusion 

In this paper, we have presented some more sophisticated applications of Lyapunov Weighted Dynamics [24] . 
In particular, we have shown that LWD can both be applied to locate atypical trajectories and to measure 
free-energy-like quantities in extended dynamical systems. This is however not yet the full-blown projects 
that one may envisage for the future, such as applications to real planetary systems, nonlinear waves, and 
fully developed turbulence. 

The method seems indeed a very powerful extension of the usual sampling of typical trajectories, 
although one must admit that it requires, compared to standard Monte Carlo algorithm, much more wisdom 
in the choice of parameters such as the number of clones, or the cloning rate. 
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Appendix A. Conserved noise 



In this section we show how to derive and implement numerically Gaussian white noises that conserve the 
total energy, and show how they lead to the microcanonical measure in the steady-state. To keep notations 
as light as possible we will choose the convention that repeated indices should be summed over. Furthermore, 
we write xf instead of ^ 



Appendix A. 1. Equations of motion 

For the most generic Hamiltonian "H, Hamilton's equations of motion read 

q^^np^ (A.l) 

Pi= - Hq, 

where we write instead of We then add a Gaussian white noise r/ and a friction that we adjust to 
keep the energy H constant. The dynamics (A.l) becomes: 

= -Hp. (A.2) 
Pi ^ ~ T-Lqi - zT-Lp^ + V^erji 
where z is such that 1-1 = 0. This implies 

H = Hq^q, + Hp^p, = Hq,Hp, + Hp^ (-Uq^ - zHp, + \/2i77,) (A.3) 
= - znl^+^/YeUp^T^, 



Finally we set 



2e ^ (A.4) 



so that the equations of motion are given by 

q, = Hp, (A.5) 



- nq^ - V2e-^Hp. + V2srj, = -Uq^ + V2e 5,k - ) Tlk 



2 ■ -Fz ■ ' (I ■ -i/i ' • — 1 -if^ n_/2 

i^pi \ i^pi 



Qik 

Let us note that g is the projector onto the energy surface in the momentum space. Indeed g satisfies 

5., - 9.k9k, - [S., - j [S,, - j _ ^S.,, - 2-^ + -J^YJ^ j(A-6) 

= ffy 

/ -"pi \ 

Furthermore, it is straightforward to check that the vector ^ = , which is normal to the energy 

\ '"■PN / 



surface in the momentum space, is in the kernel of g: 



igV)^ = g^,V, = l^S^.Hp^ - ^Ip^ j = - Hp. = (A.7) 
Hence, g is indeed the projector onto the energy surface. 

2 

When 'H{q,p) — ^ + V{q), the dynamics (A.5) becomes 

q, = p, (A.8) 

Pi = - Vg, - V2e^^^^pi + V2er]^ = - Vq^ + ^/2eg^jrij 
Pi 

where gij = Sij — ^^^. To keep notations as simple as possible, we restrict the derivation of the Fokkcr-Planck 
equation to this case, even though it extends to more general Hamiltonians. 
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Appendix A. 2. Fokker- Planck Equation and steady-state measure 

The noise in (A. 8) is multiplicative and this Langevin equation is well defined only when a prescription has 
been chosen to time-discretize g{p)i. Let us now show that under Stratonovich convention, microcanonical 
measures are steady-state of the Fokker-Planck equation. Starting from the Langevin dynamics 

i^ = h,{^,t) + V2^g,JT, (A.9) 

where Tj is a Gaussian white noise such that {T j{t)ri{t')) — 5{t — t')Sij, diffusion and drift coefficients are 
given by [63]: 



which yields for (A. 



£9kj 



da 



and 



D 



■I] — £9ik9]k 



and 



d9p. 



dpk 



^PiPj ~ ^9piPk9pjPk 
The Fokker-Planck equation then reads: 



dP 

1h 



d_ 
dq, 

d_ 
dqi 



d d_ ,d_ 



api 



■ 9pjPk 



P 



99p, 



Pk ■ 



dpj 



P 



Using the explicit expression for g^-p^ and the fact that gp^p^Pj 



dP 
~dt 



d_ 

dq,. 



dpi 



_d_ 
dp. 



d_ 
dpt 



PiP] 



'LrPr dp J 



0, we get: 
' ' P 



Let us now consider the evolution of a distribution f{q,p) that solely depends on T-L{q,p) 



Using ^ 



5/(g,p) ^ df{n) 
at dt 

weget: 
dm{q,p)] 



dq, 



p,fm- 



d_ 
dpi 



dt 



= -rm 



dn 

dq. 



-Pi 



dp/" 



' dp i 



d_ 
dpi 



d_ 
dp, 



PiPj d 



p"^ dp^ 



fin) 



dn 

dpi 



PiPj dU 
p"^ dpj 



= 



(A.fO) 

(A.ll) 
(A.12) 

(A.13) 
(A.14) 

(A.15) 

(A.f6) 

(A.17) 



We have thus shown how to construct a Langevin equation that conserves the energy and yield a uniform 
measure onto the energy surface. Let us further notes that the numerical implementation of (A. 5) could be 
rather difficult, involving a cumbersome multiplicative noise. Fortunately, there is a very simple way around 
this, exact at the order V edt but which conserves the energy exactly. We pull a random vector rj on a sphere 
of dimension N and radius VZedt, which we add to p. We then renormalise p to keep its norm constant. 
This amounts to 



P 



\P\ 



\p+^/2eAt■n\ 



(p + V2edt T] 



At first order in \/2edt, this reads 

V2edt 



P^P 



ri p 



PV 

\p\\ 



(A.18) 



(A.19) 



If one further wishes to conserve a total impulsion '^^Pi = 0, we simply replace rji by yy^ — rjj 
in the above procedure. Last, let us note that despite our manipulations, these noises remain Gaussian, as 
linear combinations of Gaussian noises. 



B The computation leading to tiie conservation of tiie energy assumes a Stratonovich convention [62]. 
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